library(smashr)
Now we consider estimating a spatially structured mean \(\mu\) = (\(\mu_1\), . . . , \(\mu_2\) ) from Poisson data:
\[y_t ∼ Pois(\mu_t), (t = 1, . . . , T)\] # Simulation 1
n = 1000
t = 1:n
spike.f = function(x) (5 * exp(-100 * round((x - 100)/10)^2) + 5 * exp(-100 * round((x - 400)/10)^2) + 5 * exp(-100 * round((x - 800)/10)^2) + 2 * exp(-100 * round((x - 500)/100)^2)
)
mu.s = spike.f(t)
# POISSON CASE
# ------------
# Scale the signal to be non-zero and to have a low average intensity.
mu.t = 1 + mu.s
# Simulate an example dataset.
set.seed(101)
X.s = rpois(n, mu.t)
# Run smash.
mu.est = smash(X.s, "poiss")
# Plot the estimated mean function (red) against the ground-truth (black).
plot(mu.t, type = "l", col="darkgreen", lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (kb)", ylim=c(0,max(X.s)))
points(X.s,col="darkgreen", cex=0.2)
lines(mu.est, col = "salmon", lwd=1.5)
Here, green is the real parameter and data. red is the esitmated mean by smash.
n = 500
t = 1:n
spike.f = function(x) ( 2 * exp(-100 * round((x - 100)/10)^2) + 2 * exp(-100 * round((x - 450)/10)^2) + 0.5 * exp(-100 * round((x - 300)/100)^2)
)
mu.s = spike.f(t)
# POISSON CASE
# ------------
# Scale the signal to be non-zero and to have a low average intensity.
mu.t = 1 + mu.s
# Simulate an example dataset.
for (i in 1:50){
set.seed(i)
X.s = rpois(n, mu.t)
# Run smash.
mu.est = smash(X.s, "poiss")
# Plot the estimated mean function (red) against the ground-truth (black).
print(i)
plot(mu.t, type = "l", col="darkgreen", lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (kb)", ylim=c(0,max(X.s)))
points(X.s,col="darkgreen", cex=0.2)
lines(mu.est, col = "salmon", lwd=1.5)}
## [1] 1
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 13
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 14
## [1] 15
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 16
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 17
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 18
## [1] 19
## [1] 20
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 21
## [1] 22
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 38
## [1] 39
## [1] 40
## [1] 41
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## [1] 42
## [1] 43
## [1] 44
## [1] 45
## [1] 46
## [1] 47
## [1] 48
## [1] 49
## [1] 50
Here, green is the real parameter and data. red is the esitmated mean by smash.
We use 20kb window size around region KCNQ2 (chr20:62,037,542).
region <- head(read.table("table.100k.window.stand.diff.cutoff.4.topwindow.across.all.mutation.type.CSV", sep=",", skip=1, stringsAsFactors = F), n=5)
dat <- read.table("table.ASCWGS_20180504.WGS1902_hg19_controls_SNV_remove_recurrent_mutation.20kwindow.csv")
for (i in 1:nrow(region)){
dat0 <- dat[dat$chr==region[i,2] & dat$start > region[i,3] - 6e6 & dat$start < region[i,4] + 6e6, ]
# Run smash.
mu.est = smash(dat0$obs, "poiss")
plot(dat0$start/1e6, dat0$obs, type = "l", col=NA, lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (Mb)", ylim=c(0,max(dat0$obs)), main = paste0(region[i,2],":",region[i,3]))
points(dat0$start/1e6, dat0$obs, col="darkgreen", cex=0.2)
lines(dat0$start/1e6, mu.est, col = "salmon", lwd=1.5)
}
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.